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ABSTRACT 

Observations show that radial metallicity gradients in disk galaxies are relatively 
shallow, if not flat, especially at large galactocentric distances and for galaxies in the 
high-redshift universe. Given that star formation and metal production are centrally 
concentrated, this requires a mechanism to redistribute metals. However, the nature of 
this mechanism is poorly understood, let alone quantified. To address this problem, we 
conduct magnetohydrodynamical simulations of a local shearing sheet of a thin, ther- 
mally unstable, gaseous disk driven by a background stellar spiral potential, including 
metals modeled as passive scalar fields. Contrary to what a simple a prescription for 
the gas disk would suggest, we find that turbulence driven by thermal instability is 
very efficient at mixing metals, regardless of the presence or absence of stellar spiral 
potentials or magnetic fields. The timescale for homogenizing randomly distributed 
metals is comparable to or less than the local orbital time in the disk. This implies that 
turbulent mixing of metals is a significant process in the history of chemical evolution 
of disk galaxies. 

Subject headings: galaxies: abundances — galaxies: ISM — galaxies: kinematics and 
dynamics — instabilities — methods: numerical — turbulence 



1. INTRODUCTION 



The spatial distribution of metals in disk galaxies is a crucial clue for understanding how 
galaxies formed and evolved over cosmic time. The past few decades have produced a wealth of 



observations of this property, including in our own Milky Way (e.g., ] 


ienry et al. 2010 Balser et 


al. 2011; Luck & Lambert 2011 Cheng et al. 


2012 


Yong et al. 2012 


'), in nearby galaxies (e.g., 


Vila-Costas & Edmunds 1992 Considere et al. 


2000 Pilyugin et al. 2004 Kennicutt et al. 2011), 


and in the high-redshift universe (e.g., Cresci et al. 


2010 Jones et al. 


2012 Queyrel et al. 2012). 



In the local universe, a variety of radial metallicity gradients in disk galaxies are seen, but they are 
generally on the order of —0.03 dex kpc -1 , with the negative sign indicating decreasing metallicity 
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at larger galactocentric radii. However, surprisingly, these gradients seem to disappear in the outer 
parts of galactic disks, where there is little star formation; moreover, the metal content is greater 
than would be expected given the amount of star formation that has taken place at these radii 



(e.g., Bresolin et al. 2009, 2012 Werk et al. 2011). High-redshift galaxies, in comparison, are far 
less regular. Their metallicity gradients range from negative ones significantly steeper than those 
found locally, to completely flat or even positive. Any successful theory of galactic evolution must 
be able to reproduce these observations. 

Several processes play an important role in regulating the spatial variation of metals in disk 
galaxies, and these have been demonstrated by either chemical evolution models or hydrodynamical 
simulations. The enrichment of metals in the interstellar medium (ISM) is dominated by star 
formation and subsequent stellar mass loss, and thus depends on the star formation law (e.g., 



Phillipps & Edmunds||1991 


). Metals are diluted by infalling gas from outside galaxies (e.g., Tinsley 


& Larson||1978 |Chiosi||198 


D Matteucci & Frangois 1989; Chiappini et al.||1997; Prantzos & Boissier 


2000 Chiappini et al.|2001 


) . Radial inflow of the gas within the disk of a galaxy redistributes metals 


(Mayor & Vigroux 1981 L 


;acey & Fal 


1 1985, Pitts & Tayler 1989 


Gotz & Kbppen 1992 Portinari 


& Chiosi 2000 Spitoni & Matteucci 


2011 Bilitewski & Schonricr 


l 2012); galaxy interactions are 



especially effective in inducing large-scale inflow and flattening metallicity gradients (Rupke et al. 



2010 Perez et al. 2011; Torrey et al., in preparation). Turbulence associated with the viscous 



evolution of gas disks can also redistribute metals (Clarke 1989; Sommer-Larsen & Yoshii 1989 
1990; Tsujimoto et al.|1995 Thon &: Meusinger|1998 ). Beyond these gas-dynamical processes, radial 
migration of stars can alter stellar metallicity distributions independently of processes affecting the 
gas phase ( |Sellwood fc Binney|[2002| |Roskar et aL]|2008a|bt |Schonrich Binney|[2009l ) . 



However, the strength and relative importance of these processes remains very poorly under- 
stood. Semi-analytic chemical evolution models generally parameterize each process and then tune 
the parameters in an attempt to provide an acceptable match to observations. However, the large 
numbers of parameters involved means that even a good fit to the data may not be unique, and 
the need for fine-tuning means these models have limited predictive power. Moreover, the param- 
eterizations used in the models may not be accurate. For example, turbulent mixing is usually 



treated by adopting an a prescription for the turbulent transport of angular momentum (Shakura 
Sunyaev| |l~973), and assuming that the transport coefficient for metals is the same. Under these 
assumptions turbulent mixing is unimportant unless a is so large that the viscous diffusion time 
becomes comparable to the gas depletion time. However, neither the assumption that turbulent 
transport of metals can be approximated with an a prescription, nor that a for this process is 
the same as that for the angular momentum, are physically well-motivated. Numerical simula- 
tions of metal transport that evolve galaxies over cosmological times are unfortunately little better, 
because their limited resolution means that they must also adopt parameterized treatments of un- 
resolved processes. For example, most smoothed particle hydrodynamics (SPH) simulations allow 



no chemical mixing at all between SPH particles (Wadsley et al. 2008), or at best treat mixing 



approximately using a parameterized subgrid recipe (Shen et al. 2010). Eulerian simulations, in 
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contrast, dramatically over- mix when their resolution is low. 



The approach we take in this paper is quite different, and complementary to chemical models 
and large-scale cosmological simulations. We isolate a single process: turbulent mixing within a 
galactic disk. Our goal is to provide a first-principles calculation of this process, which can in turn 
provide a physically-motivated, parameter-free prescription that can be used in chemical evolution 
models or lower resolution simulations. In order to achieve this goal, we simulate turbulent mixing 
in a portion of a galaxy at very high resolution, including physical processes that are too small-scale 
to be resolved in cosmological simulations, and we perform a resolution study to ensure that our 
results are converged. The previous work that most closely matches ours in philosophy and overall 



approach is that of Mac Low & Ferrara (1999) and Fragile et al. (2004), who used high resolution 



simulations of isolated portions of galaxies to study mixing of supernova ejecta with the ISM as 
galactic winds are launched. Here we perform a similar calculation for turbulent mixing within 
disks. 



There exist many sources that can drive turbulence in the ISM (see Elmegreen & Scalo 2004 
and Scalo Elmegreen|2004 and references therein). In earlier work, de Avillez & Mac Low (2002) 
studied the properties of turbulent mixing driven by supernova explosions. Here, we instead fo- 



cus on turbulence driven by thermal instability (Field 1965; Field et al. 1969). Our motivation is 



two-fold. First, the flat metallicity gradients seen in outer disks presumably call for some sort of 
mixing process to operate at large galactic radii, where star formation is limited and thus supernova 
explosions are extremely rare. Second, even in places where supernovae do occur, thermal insta- 
bility is also present and will drive turbulence; indeed, thermal instability is essentially inevitable 
anywhere the ISM is dominated by atomic hydrogen, which is the case for most galaxies over the 
great majority of cosmic time. Supernovae will only enhance the turbulence compared to what we 
find, and thus our results should be viewed as a minimum estimate of the turbulent mixing rate. 

In the sections that follow, we describe our numerical method and present the results of the 
simulations. We then quantify the results in a form appropriate for use in chemical evolution models 
and discuss their implications. Finally, we summarize the results and conclude. 



2. NUMERICAL MODELING 

To study turbulent mixing, we adopt a thin-disk, local-shearing-sheet model similar to |Kim 



Ostriker| d2002~l ) (see also |Kim & Ostriker|[2006l ) and equip it with approximate heating and cooling 



processes (see also Kim et al. 2008, 2010) and metal tracers. We also investigate the effects of 
spiral shocks and magnetic fields on turbulent mixing in our models. In the following sections, we 
describe our simulation methodology and setup in detail. 
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2.1. Governing Equations 

2.1.1. Magnetohydrodynamics 



Using the local-shearing-sheet approximation (Goldreich & Lynden-Bell 1965), we consider a 
small region distant from the center of a vertically thin disk galaxy. This region is centered at 
the potential well of a background stellar spiral arm and co-rotating with the arm at its angular 
pattern speed Q p . One can define a coordinate system for this region by x' = R — Ro and y' = Rcf), 
where (R,4>) are the polar coordinates rotating with the spiral arm and (R,<j>) = (Ro,0) is the 
location of the point of interest along the stellar spiral arm. The velocity field of the gas flow in the 
rotating frame is denoted by v. Since the differential rotation profile of a disk galaxy Q = 
is usually a known or given function, we are more interested in the relative velocity of the gas 
u = v — v c with respect to the circular velocity v c = R(£l — Qp)e y i in the rotating frame. By 
assuming x' <C Ro, y' <C Ro, u x > -C -Ro^O) and u y > <C -Ro^o> where flrj = ^{Ro), and expanding 
the magnetohydrodynamical equations to first order in x' , y' , u x i, and u y >, the continuity, the 
momentum, the energy, and the induction equations become 



dt 

d_ 
di 



as +v VS + SV-u = 0, (1) 



^ + v • Vu = qon u x ,e y , -2ft xu- ±Vp - V ($ s + $ g ) + ^ J x (B + B) , (2) 



8A 

dt 



|f + v Ve + eV-u = -pV-u + H, (3) 
+ v c • VA = qi QoA y/ e x , + u x (B + B) , (4) 



respectively. The primitive variables for which we solve the above equations are the gas surface 
density S, the gas relative velocity u as defined above, the thermal energy density of the gas e 
(i.e., internal energy per unit surface area), and the magnetic vector potential A. The magnetic 
field B is then calculated by B = V x A. The remaining quantities in the above equations are the 
dimensionless shear parameters 



RdQ 
TldR 



(5) 

R=Ro 



1 dfl (fi - Op 



n dR 



= 90-1 + ^, (6) 

R=Ro U 



the gas pressure p, the gravitational potentials due to the stellar spiral arm and the gas itself, <3? s 
and respectively, the electric current density J = (V x B)//io, where fio is the permeability, and 
the net heating rate per unit surface area H. We impose a constant external azimuthal magnetic 
field Bo = Boe y i when we consider a magnetized diskQ 

Equations @-@ are written in vectorial forms and thus are readily transformed into different 
coordinate systems. One particular choice is to rotate the (x',y r ) system counterclockwise by the 



1 The induction equation (4 1 requires that Bo || v c ; otherwise, an additional term v c x Bo should be included. 
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pitch angle i into the (x,y) system such that the new a;- axis and y-axis are perpendicular and 



parallel to the stellar spiral arm at the origin, respectively (see Figure 1 of Kim & Ostriker 2002). 
We employ the tightly-wound approximation so that sin i ~ i <C 1 can be deemed a small quantity. 
In these considerations, Equations and ^ remain unchanged while Equations ^ and Q 
become, by also retaining sini to only first order, 



h + v • Vu = q n u x e y - 2fi x u - ^Vp - V ($ s + $ g ) + ± J x (B + B) , (7) 
' ' A + v c • VA = gift [{Ax sin i + A y ) e x - A y e y sin i] + u x (B + B) , (8) 



at 



respectively. The circular velocity and the external magnetic field in this tilted frame can be 
approximated by 

v c « v e x sin i + (v - qiSl %) e y , (9) 

B » -B e x sin « + B e y , (10) 

respectively, where = i?o (f^o — ^p)- 



2.1.2. Forcing Driven by the Stellar Spiral Arm 



The advantage of aligning our coordinate system with the background stellar spiral arm is that 
the gravitational potential of the arm in this system can be approximated as periodic in x while 
weakly varying in y, to first order flRoberts|[l969"t |Shu et al.|l973| |Kim & Ostrike71|2002[ ) : 



$ s {x) ~ $ocos 



2irx 



where 3>o < is a constant, 



L 



2ttRq sin i 



in 



(11) 



(12) 



is the radial spacing between adjacent spiral arms, and m is the multiplicity of the arms. The 
strength of the spiral forcing is measured in terms of the local centrifugal acceleration: 



F 



m 



sin i 



(13) 



Note that according to Equation (12), the local-shearing-sheet approximation L x <C Rq requires 
that sini/m <C 1, which is automatically satisfied by the tightly-wound approximation sini <C 1. 



2.1.3. Self-gravity of the Gas 



In this work, we are only interested in large-scale mixing of metals and ignore self-gravity of the 
gas. It will be required, though, in a subsequent paper where we investigate the metal abundances 
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in precursors of molecular clouds. Therefore, we list the equation governing self- gravity of the gas 
in this section for completeness. 

The quantity $ 9 in Equations ^ and Q represents the gravitational potential of the gas. 
The potential is calculated by solving the Poisson equation for a razor-thin disk 

V 2 $ 9 = 4ttG£5(z), (14) 

where G is the gravitational constant and 5{z) is the Dirac delta function. We discuss the corre- 



sponding numerical method for its solution in Section 2.3 



2.1.4- Thermodynamics 

To understand the effects of thermal instability on the mixing of metals, we compare thermally 
stable disks with thermally unstable ones. For the former, we use the isothermal equation of state 
p = c^S, where c s is the isothermal speed of sound, and in this case, the energy equation Q 
is not required and we do not solve it. For the latter, we adopt the adiabatic equation of state 
p = (7 — l)e, where 7 is the two-dimensional adiabatic index, and include the heating and cooling 
of the gas such that the disk is thermally unstable. 

For a prescription of the heating and cooling rates, we start from the approximate functions 
suggested by Koyama & Inutsuka ( 2002[ )p| The net rate of heat loss per unit volume is pC = 
n 2 A — riT, where n is the number density of gas particles and 



r 

A(T) 



2.0 x 10 



10 7 exp 



-26 



erg s , 
1.184 x 10 5 



+ 1.4 x 10" 



Texp 



92 



cm 



(15) 
(16) 



I V T + 1000 

in which T is the temperature in Kelvins. To obtain the heating rate per unit surface area H, we 
need to integrate pC in the vertical direction, and the vertical structure of the gas is required. For 
simplicity, we assume that the gas is vertically isothermal and the number density is approximated 
by n(z) ~ no exp [—z 2 /2H 2 ^ , where no is the number density in the mid-plane and H is the vertical 
scale height. We further assume that the vertical motion of the gas is dominated by non-thermal 
processes, at least when the gas temperature is low, such that H is constant, instead of depending 
on T. Therefore, 



n 



pC dz 



pm u 



1 



1 



pm u 



A(T) 



(17) 



where we have used j ndz = E/um u , and u and m u are the mean molecular weight and the atomic 
mass, respectively. To complete the system, we use the ideal-gas law p = Efc^T '/ ' pm u for the 
temperature T, where ks is the Boltzmann constant. 



2 The cooling function published in 



Koyama & Inutsuka 



( 2002 1 contains two typographical errors and has been 



corrected by Nagashima et al. ( 2006 1 
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2.1.5. Metal Tracers 

Finally, we assume the heavy metals in the disk follow the velocity field of the gas component 
and model them as tracer fluids. Therefore, the surface density of any given metal T,x satisfies 

-^ + v-VS x + S x V-u = 0. (18) 

The concentration of the metal c with respect to the total local mass can then be derived by 
c = S x /S. 



2.2. Initial and Boundary Conditions 

The computational domain we consider is a square sheet of size L, the radial spacing between 
adjacent spiral arms defined in Section 2.1.2 In the following subsections, we discuss the initial 
and boundary conditions we adopt for the gas and the metals. 



2.2.1. The Gas and the Equilibrium State 

We initially set the gas to be uniform (E = So), isothermal (T = To), and moving along 
with the galactic circular motion (u = 0) such that it is at an equilibrium state when the spiral 
forcing is not present (<3?o = F = 0). Our adopted initial density So can be more physically 
motivated by considering the corresponding Toomre Q parameter for the gas Qq = «;c S) q/7tGSo, 
where k is the epicycle frequency of the disk at R = Rq and c S) o is the initial speed of sound. The 
epicycle frequency k as well as the shear parameters qo and q\ defined in Equations ^ and ([6| are 



determined by the rotation profile Q(R). As in Kim & Ostriker (2002), we assume a flat rotation 



curve (RQ = constant) near R = Rq and a pattern speed of Vt p = Qq/2 and thus k = v2^o> Qo = 1> 
and q\ = 1/2. In physical units, the initial surface density is then 

S = (19 M pc- 2 ) Q 1 ( ^ T ) ( P^—r) . (19) 

V OP V26kms- 1 kpc-V Vr.Okms- 1 ; V ' 

On top of the equilibrium state, we perturb the velocity field by a white noise of magnitude 10 _3 c Sj o 
to seed the instabilities, if any, of the system. 

When we consider the isothermal equation of state, our initial isothermality of the gas is 
automatically guaranteed and preserved. In this case, only the constant speed of sound c s = c Sj o 
needs to specified. On the other hand, when we consider a non-isothermal disk with heating and 
cooling processes, an initial thermal equilibrium of the gas is also required. This is equivalent to 
setting % = at S = Sq and T = Tq, i.e., zero net heating rate at the initial state. Given a 



surface density S, Equation (17) can be used to solve for the corresponding temperature T such 



that % = 0. With the values of the physical parameters considered in this work (see Section 2.2.3 
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Fig. 1. — Pressure-density diagram for non- isothermal gas in our models. The solid line indicates 
the states at thermal equilibrium, the dotted lines show the isotherms, and the dashed lines denote 
constant Toomre stability parameter for the gas. 



and Table [T]) , Figure [T] plots the curve for the states at thermal equilibrium in a pressure-density 
diagram; notice the region in S/ fim 

10 2i_ 1Q 22 

cm where the slope of the curve is inverted, the 



classical condition for thermal instability to occur (Field 1965). Therefore, a given value of initial 
Toomre stability parameter for the gas Qo uniquely specifies the initial state of the gas (Eo,Tb). 

For magnetized disks, we impose initial uniform azimuthal fields through the gas by setting 
Bo / and A = throughout. The strength of these imposed fields can be gauged by the 
corresponding plasma beta parameter, (3q = 2/ioPo/ Bq, which is the ratio of the initial thermal 
pressure po to the initial magnetic pressure. 

Since our system is driven by a forcing periodic in the x direction with a wavelength of the 



inter-arm spacing L (Equation ( 11 )), it is expected that the response of the system also be periodic 



in x with the same wavelength. The system is also sheared in the y direction, though, which is 
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manifested by the term — qi^Qxe y in Equation (j9j). So strictly speaking, the boundary conditions in 
the x direction for the system should be sheared periodic, which can be expressed mathematically 
in the form f(x + L, y) = f(x, y + qiQ^Lt), where f(x, y) is any dynamical field in question (Hawley 



et al.| 1995 Brandenburg et al. 1995 |Kim &: Ostriker 2002), except the metal tracer fields (see 
below). As for the y direction, we adopt normal periodic boundary conditions f(x, y + L) = f(x, y) 
for convenience. 



2.2.2. The Metal Tracers 

Given that we model the metals as passive scalar fields, they effectively act like dye in a 
flow. We can in fact inject metals anywhere, anytime, and at any rate, to study their diffusion 
process. Since most star formation occurs along spiral arms, constantly producing metals that drift 
downstream towards the next spiral arm, we are most interested in the mixing of metals within 
one passage between adjacent arms. This motivates us to employ inflow boundary conditions at 
the left (x = —L/2) and outflow boundary conditions at the right (x = +L/2) for the metal tracer 
field "Ex- The boundary conditions in the y direction remains periodic. 

The spatial distribution of the newly produced metals from the previous generation of star for- 
mation may be arbitrary. To be as general as possible, we constantly inject a sinusoidal distribution 
of unit amplitude from the left: = —L/2,y) = sin(2-7ry/Ai n j), where Ai n j is the wavelength 

of the distribution. Once the wavelength dependence of the mixing process is deciphered, an ar- 
bitrary distribution of metals can be analyzed by Fourier decomposition. In all of our models, we 
simultaneously evolve four species of metal tracers with A; n j = L, L/2, L/4, and L/8, for which we 
denote X by 1, 2, 3, and 4, respectively. 

At the equilibrium state when there exists no driving force in the system, the gas flows az- 
imuthally at circular velocity, and so do the metals. Given that the azimuthal direction in our 
computational domain is tilted from the y-axis by the pitch angle i, we set the initial condition for 
the metal tracers as 

2vr / x + L/2 s 



E x (x,y) = sin 



(20) 



2.2.3. Physical Parameters and the Models 



The values of the physical parameters we adopt and keep constant across different models are 
listed in Table [TJ Most of these values match those used by Kim & Ostriker ( 2002 ) for comparison 
purposes. The additional two parameters, mean molecular weight \i and vertical scale height H, 
come into the system only via the thermodynamics. For simplicity, we set \i = 1 throughout]^] The 



3 A more realistic value should be fi ~ 1.3. But this difference does not significantly change our results. 
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vertical scale height is estimated by noting that H = c StZ /v in a vertically isothermal gas with 
effective vertical speed of sound c S)Z and vertical frequency v\ in the Solar neighborhood, v ~ 2Q 

t 7.0 km s^ 1 , we calculate H 



(Binney & Tremaine 2008). By assuming c 



95 pc. 



We conduct separate simulations with each combination of three physical effects — spiral 
forcing, thermal instability, and magnetic fields — to study their influence on the gas dynamics and 
more importantly the mixing of metals. The eight resulting models and their numerical resolutions 
are listed in Table[2j We show in Appendix[B]that this resolution is sufficient to achieve numerically- 
converged results. If a model includes spiral forcing, we gradually increase its strength for a duration 
of 10P to F = 3%, after which the strength remains constant. When considering the isothermal 
equation of state, we use a speed of sound of c Sj o = 7.0 km s _1 and thus an initial gas surface 
density of So = 13 Mq pc~ 2 . When considering thermally unstable gas, we adopt a two-dimensional 
adiabatic index of 7 = 1.8, which is taken to be the limiting value of a strongly self-gravitating 



disk of monatomic gas (Gammie 2001). The initial thermal equilibrium in this case requires that 
S = 12 Mq pc' 2 and c Sj0 = 6.4 km s" 1 . For magnetized disks, we set the initial plasma beta to 
be A) = 2. 



2.3. The Pencil Code 



We use the Pencil Cod^] to solve our system of equations discussed above. It is a cache- 
efficient, parallelized code optimal for simulating compressible turbulent flows. It solves the MHD 
equations, among others, by sixth-order finite differences in space and third-order Runge-Kutta 



steps in time, attaining high fidelity at high spectral frequencies (Brandenburg 2003). Although 
the scheme is not written in conservative form, conserved quantities are monitored to assess the 
quality of the solution. 

Several diffusive operations are employed in order to stabilize the scheme. We use hyper- 
diffusion in all the four dynamical Equations ([!]), ^ , ([7]) , and Q to damp noise near the Nyquist 



frequency while preserving power on most of the larger scales (Haugen Sz Brandenburg 2004; Jo- 



hansen k Klahr|2005 


). Shocks 


et al. 


2004 


Lyra et al. 


2008 


)• 



maintain roughly the same strength of diffusion at the grid scale (see Appendix [A]) . Finally, all 
the advection terms of the form (vo + u) • VQ, where vo = voe x sini + voe y (see Equation Q) 
and Q is any state variable, are treated by fifth-order upwinding to avoid spurious oscillations near 



stagnation points (Dobler et al. 2006). 



Since a local shearing sheet is considered, we need to handle the sheared advection, the bound- 
ary conditions, and the Poisson equation with care. The sheared advection terms of the form 



-qi^loxdyQ are directly integrated by Fourier interpolations (Johansen et al. 2009) in order to 



4 The Pencil Code is publicly available at http://code.google.eom/p/pencil-code/. 
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relieve the time step constraint from shearing velocity (Gammie 2001) and eliminate the artificial 



radial dependence of numerical diffusion (Johnson et al. 2008). The sheared periodic boundary 
conditions discussed in Section 2.2 are similarly implemented with Fourier interpolations. The 



Poisson Equation (14) is solved by fast Fourier transforms in sheared Fourier space in which the 



fields are strictly periodic (Johansen et al. 2007) 



We have implemented the approximate net heating function (Equation (17)) in the Pencil 
Code. Since the thermal and the dynamical timescales can be quite different, we operator split 
this term in the energy Equation ([3]). Because this heating and cooling process only depends on 
local properties, the resulting differential equation is ordinary and can be integrated independently 
at each cell. For these integrations, we adopt the fifth-order embedded Runge-Kutta method with 



adaptive time steps (Press et al. 1992 ). Since most computational cells are near thermal equilibrium, 



they require only one or two iterations to match the hydrodynamical time step. Integration of the 
remaining few cells that have shorter thermal times has negligible computational cost. 

With our highest resolution of ~1.5 pc per cell, we are still not able to resolve the thermally 
stable cold phase of the gas (see Figure [T]). Therefore, the densest cells tend to overcool and lose 
pressure support to their surroundings. To ensure the Jeans length is properly resolved and avoid 



artificial fragmentation (Truelove et al. 1997) later when we include self-gravity of the gas, we 



impose a floor to thermal energy density e in accord with the local surface density £ such that the 
condition of at least four cells per Jeans length is satisfied: e > AGT?h/^{^ — 1), where h is the cell 
size. This is in effect a modification to the cooling function at low temperatures. Given that our cell 
size is marginally close to resolving the stable cold phase and our main purpose is to demonstrate 
if thermal instability can drive effective chemical mixing, we omit any further consideration of 
sub-scale physics and leave this compromise as a caveat. 



3. STATISTICALLY STEADY STATE 

All of our models attain statistically steady state within a few local orbital periods. In this 
section, we report the state of our simulations at this stage. 



3.1. Gas Dynamics 

Figures [2] and [3] show the snapshots of the density field for the non-magnetic and magnetic 
models, respectively. For models Control and M, since there exists no driving force in the system 
(i.e., neither spiral forcing nor thermal instability), no interesting feature occurs and these models 
serve as control simulations; the initial perturbation propagates as sonic waves and remains small 
in amplitude. For model F, isothermal disk with spiral forcing, a spiral shock with little azimuthal 
variation forms and locates slightly upstream of the potential well of the forcing. As can be more 
clearly seen in the y-averages of the gas properties plotted in Figure [4j this resembles the classical 
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solution of iRobertsI (|1969l) 



When the disk is thermally unstable, the gas spontaneously breaks into two phases, a warm 
phase of diffuse gas with high volume filling factor and a cold phase of dense gas with filamentary 
structures, as evident in models T and TF shown in Figure [2] Model T is simply the standard 
two-phase model regulated by thermal instability ( Field et al,| [T969). The two-phase medium is 



turbulent but statistically steady, in which the two phases are in approximate pressure equilibrium. 
As demonstrated by model TF, although the spiral forcing tends to concentrate material into the 
potential well, the spiral shock seen in isothermal disks is suppressed by the turbulent medium 
(Figure Q . 

The existence of magnetic fields creates an interesting structure in the gas. For model MF, 
magnetized isothermal disk with spiral forcing, two discontinuities parallel to the spiral arm occur as 
shown in Figures [3] and |4j The gas remains at roughly the initial state in between the discontinuities, 
where the magnetic field lines are compressed and magnetic pressure is increased. While the left 
discontinuity is quite stable, the right discontinuity becomes wobbly after the spiral forcing reaches 
its maximum strength. Waves are produced in the process and they propagate throughout the 
disk. This behavior, however, does not continue to develop in magnitude and drive the gas into a 
turbulent state. 

Finally, the magnetized, two-phase, turbulent medium is rather similar to its non-magnetized 
counterpart, as can be seen by comparing model MT shown in Figure [3] with model T shown in 
Figure [2] As in the non-magnetized case, the presence of thermal instability significantly weakens 
the spiral shock, as demonstrated by model MTF shown in Figures [3] and |4j 



3.2. Metal Tracers 

With a statistically steady state of the gas established for each model listed in Table [2j we 
turn to observe how metals would be transported in each flow. Figures [5}{8] show snapshots of 
the metal tracer fields with an injection wavelength of A; n j = L or L/2 for each model at time 
t = 15P. Since models Control and M remain at their initial equilibrium states, as described in 



Section 3.1, we expect the metal tracers do the same. In this case, metals should move in the 
azimuthal direction (which is tilted to the right at an angle i with respect to the y-axis) without 
any noticeable diffusion^] and our simulations have passed this benchmark as demonstrated in the 
snapshots. 

When the spiral forcing is present in our non-magnetized isothermal disk (model F), the tracer 
field is deformed in the x direction, which is perpendicular to the spiral arm. In accord with the 
velocity field shown in FigureHl the metal distribution is first rarified and its wavelength is increased, 



We note that molecular diffusion is too small to drive metal diffusion on galactic scale, and this process is 
obviously ignored in our models. 
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Fig. 2. — Snapshot of the density field for each non-magnetic model at t = 15P. Isothermal disks 
are in the left column, while thermally unstable disks are in the right column. The top row has no 
spiral forcing, while the bottom row does. The color scales are set the same for all the panels. 
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Fig. 4. — The y-averages of the density, velocity, and magnetic fields as a function of x. These 
profiles are also time averaged from t = 10P to t = 15P. 
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Fig. 7. — Snapshot of the metal tracer field with an injection wavelength of Ai n j = L/2 for each 
non-magnetic model at t = 15P. The arrangement is the same as in Figure [2] 




Fig. 8. — Snapshot of the metal tracer field with an injection wavelength of Ai n j = L/2 for each 
magnetic model at t = 15P. The arrangement is the same as in Figure |3| 
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as the gas approaches the spiral shock. As the gas passes through the shock, the metal tracer 
amplitude is increased to ~5.5 and the wavelength is reduced. Finally the metal field is rarified 
again to regain the original distribution towards the right boundary. A similar process occurs in our 
magnetized isothermal disk (model MF), except that in this case there exist two discontinuities and 
propagating waves generated by the wobbly shock on the right. The waves, however, do not have 
enough strength to stir the metals, and the metals again tend to retain their original distribution 
after crossing the right shock. Therefore, the spiral forcing in our isothermal disks, either non- 
magnetized or magnetized, does not have a noticeable net effect on metal distribution after one 
passage of the spiral arm. 

A completely different scenario for transporting metals occurs in thermally unstable disks. 
As evident in model T, the turbulence driven by thermal instability significantly churns up the 
metals, and within only a few wavelengths in distance, the original sinusoidal distribution cannot 
be discerned anymore. To quantify this process, we compute the power spectra of the metal tracer 
fields at our final times, 

P x (k) = \t x (k)\ 2 , (21) 

where £x is the Fourier transform of a given metal tracer field. We compute the Fourier transform 
and thus the power spectrum only for gas in the downstream region, defined as the region x > 
for the runs without spiral arm forcing, and as the region beyond the spiral shock for runs with 
forcing. For models Control and M, the power spectrum is simply a 5 function at the injection 
wavelength, while for model F it is a 5 function at a wavelength smaller than the injection scale 
(due to compression of the wavelength in the spiral shock). Model MF is not quite a 5 function, 
but is nearly one. In contrast, Figure [9] shows the results for models T, TF, MT, and MTF. We see 
that the initial large-scale variation in metal density is redistributed to many different scales by the 
turbulence, and the resulting distribution becomes in fact white noise. Furthermore, this process 
does not depend on the injection wavelength, at least in the range L/8 < Ai n j < L simulated in our 
models. 

The mixing of metals driven by thermal instability is equally effective among all of our ther- 
mally unstable disks. By comparing Figure [6] with Figure [5] (or Figure [8] with Figure [7]), the metal 
tracer fields do not exhibit noticeable differences between models MT and T, indicating that mag- 
netic fields play little role in limiting the redistribution of metals by the turbulence. As shown in 
the same figures, although the mixing of metals is less effective in the pre-shock region when spiral 
forcing is present, the mixing process is significantly accelerated near the shock front, resulting 
again in white noise in the aftershock region (Figure [9]). Therefore, we determine the turbulence 
induced by thermal instability in our models is the only major mechanism in driving mixing of 
metals, and this mechanism can effectively redistribute metals into white noise within less than 
inter- arm distances. 
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Fig. 9. — Power spectra of the metal tracer fields with different injection wavelengths Ai n j in the 
downstream region of our thermally unstable disks. The downstream region is selected as the 
domain x > if no spiral forcing is present, or the aftershock region determined from Figure [4] if 
spiral forcing is present. The spectra are smoothed in 20 radial logarithmic bins and normalized 

2 

by J Sx(^) d/c = 1, where Sx(fc) is the Fourier amplitude of the tracer field £x at wavenumber 
k. The values of k at which the power spectra begin to decrease are the Nyquist frequencies in the 
simulations; the turndown is at lower k in model MTF due to the lower resolution of this model. 
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4. QUANTIFYING THE MIXING PROCESS 

4.1. Diffusion Coefficients 

Having established that the turbulence driven by thermal instability is the primary mechanism 
for mixing the metals, irrespective of the existence of spiral forcing and/or magnetic fields, we 
focus our attention on our model T and attempt to quantify the mixing process. It is not clear 
yet if this turbulent mixing can be described as a diffusion process and, if so, what diffusion 
coefficient describes it. In principle, these questions could be investigated by, for instance, the 



recently developed test-field method (Brandenburg et al. 2009 Madarassy &; Brandenburg 2010) 



However, we defer this more comprehensive analysis and present a toy model for an order-of- 
magnitude estimate of the mixing strength and timescale. 

We start by considering an observer who co-moves with the background advection (Equa- 
tion Q) and measures the distribution of the metals in the y direction, after the turbulent flow 
has reached its statistically steady state. We define t = (x — Xq) /v c>x as the advection time, where 
xq = —L/2 is the x coordinate of the left boundary, and at every position x and thus advection 
time t we compute the one-dimensional Fourier transform Tj% of Ex in the y direction. From this 
we compute the one-dimensional power spectrum 

Pl(t,k y ) = \t\{i,k y )\\ (22) 



We plot the result at several values of t for Si in Figure 10 At small t the metal distribution 
is very close to the sinusoidal one injected from the left boundary, and thus the power spectrum 
shows that almost all the power is in a single, long-wavelength mode. As the observer moves to 
the right, turbulent mixing redistributes the metals into many different scales while attenuating 
the amplitude of the initial distribution in the process. In time, the distribution becomes white 
noise and the power at all scales decays roughly synchronously while the gas flows towards the right 
boundary. 

If the process of redistributing metals were truly a diffusion process in the y direction of the 
observer's frame, then the distribution of metal tracers as a function of t and y would obey 

8L X B/ D ^ 



dt dy \ By J 



where D is the diffusion coefficient j^] If we assume D is a constant, the solution to Equation (23) 
with an initial sinusoidal distribution of wavenumber /c; n j = 27r/Ai n j is 

T. x {t,y) = V>(i)sm(/c inj y), (24) 



Note that we have not shown that turbulent mixing of metals really is a diffusion process, and indeed it is 
probably more complex than that. However, parameterizing in terms of a diffusion coefficient still provides a useful 
guide to the strength of the effect. 
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Fig. 10. — Power spectrum in the y direction at different advection times i, of the metal tracer 
field with an injection wavelength of Ai n j = L from model T. Each spectrum is averaged over 10 
snapshots at regular (physical) time interval from t = 6P to 15P. 
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where 

^(t)=^exp(-t/r D ), (25) 

in which ipQ = 1 is the amplitude of the injected distribution from the left boundary and td = 
1/Dfe^ n j is the time constant. Therefore, the power of the metal distribution at the injected wave- 
length decays exponentially according to ff-(t, Mnj) = ^ 2 {t) = V'o ex P ( — 2t/T£>). 



Figure 11 plots the power of the metal distribution in the y direction at the injection wavenum- 
ber Px(t> kinj) as a function of the advection time t in our model T. Two distinct stages of expo- 
nential decay can be seen for each metal tracer field, the first of which is steeper than the second. 
The transition time to between the two stages marks the time required for the metal distribution 
to become white noise, i.e., well mixed due to the turbulence. The shorter the wavelength of the 
injected distribution Ai n j, the faster the metals are mixed. With to identified for each tracer field, 
the decay of the power at each stage can be fitted separately by an exponential function as shown 



by the straight lines in Figure 1 1 , and the resulting slopes can be converted into the decay time 
constant tjj and the diffusion coefficient D by the formulae given above. The numerical values of 
to, td, and D for each Ai n j in our model T are listed in Table [31 



4.2. Implications for Chemical Evolution of Disk Galaxies 

The timescales we have measured in our model T indicate that turbulent mixing of metals 
driven by thermal instability is an efficient process, especially for the first stage discussed above. 
The time required to eradicate kpc-scale variations in metals, i.e., to in Table [3| is short compared 



to the orbital timescale P. As we have mentioned in Section 2.2.2, the sinusoidal distribution of 
metal tracers we inject along the left boundary can be considered as the metal enrichment powered 
by supernovae along a spiral arm. If the characteristic separation between star forming sites along 
a spiral arm is on the order of one kpc, the metals they produce will become well mixed after 
~30 Myr of advection. 

The second stage of turbulent mixing we see in our models is probably more relevant to long- 
term chemical evolution of disk galaxies. At this stage, metals are randomly distributed in the 
ISM and are constantly transported by large-scale convective motions of the gas. Like what occurs 
at the first stage, the signals of the metal variations at all wavelengths decay exponentially with 
time, although somewhat more slowly. The decay time constant td is on the order of ~100 Myr 
and is relatively insensitive to wavelength. Therefore, if there exists any metallicity gradient on a 
kilo-parsec scale, the gradient should be e-folded in roughly the same timescale, and this timescale 
is comparable to but still less than the orbital timescale of the galaxy. 

In this regard, turbulent mixing of metals should be an important physical process in chemical 
evolution of disk galaxies, and should be included in chemical evolution models. Although the toy 
model for turbulent mixing we presented in the previous section may not quantitatively describe the 
full dynamics of metal transport in turbulent ISM, we should have captured an order-of-magnitude 
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Fig. 11. — Power of the metal tracer field at the injection wavelength Ai n j as a function of the 
advection time t in model T. Four tracer fields with different A; n j are shown by solid lines. The 
power is averaged over 10 snapshots at regular (physical) time interval from t = 6-P to 15P. The 
dotted lines are the regression fit of an exponential function for the first stage of the mixing process, 
while the dashed lines are that for the second stage of the process. 
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estimate of the mixing strength. The diffusion operation along with the diffusion coefficient we 
have measured may serve as a simple starting point for a sub-grid recipe in chemical evolution 
models and cosmological simulations. 

We note that the diffusion coefficient of the second stage we find for kpc-scale distributions 
is on the same order of c s fiH ~ 0.7 kpc 2 Gyr -1 , even though the gas disk is not self-gravitating 
and presumably has a low Shakura Sz Sunyaev (1973) a parameter. In fact, a = {Lu x u y ) /Tiqc^ 
in our model T is about 10~ 2 , where () denotes the spatial average of the quantity enclosed. This 
demonstrates that the transport of metals does not strictly follow the viscous evolution of the gas 
disk. The convective motion of the gas can actually carry the metals over larger distances than a 
pure viscous stress allows. Therefore, the assumption of the same a prescription for both the gas 
and the metals in a chemical evolution model is not correct. 



5. CONCLUSIONS 

In this work, we simulate a local patch of a vertically thin disk galaxy and study the transport 
of metals with a variety of physical conditions. Specifically, we investigate the ability of thermal 
instability, spiral shocks, and/or magnetic fields to homogenize metals. We find that turbulence 
driven by thermal instability is especially effective in mixing the metals, regardless of the presence 
or absence of spiral shocks and magnetic fields. 

We observe two different modes of turbulent mixing in our thermally unstable disks. The first 
mode is for the turbulent gas to stir large-scale variations of metals into a random distribution. The 
timescale for this mode is short compared to the local orbital time in the galaxy, and this mode may 
contribute to obliterate the chemical inhomogeneities introduced by star forming activities along 
spiral arms. The second mode is for randomly-distributed metals to be continually homogenized 
over time by the turbulence. We find the timescale for this process is relatively insensitive to 
wavelength and is on the order of half the orbital timescale. This mode of turbulent mixing, 
therefore, should be of significance in reducing the metallicity gradient in a disk galaxy. 

We find that turbulent mixing of metals driven by thermal instability is more efficient than 



what a simple Shakura & Sunyaev (1973) a prescription of viscosity for the gas would suggest. The 
convective motion of the turbulent gas can in fact transport metals over larger distances, especially 
for kpc-scale variations. The dynamics is perhaps more complicated than ordinary diffusive trans- 
port with a constant coefficient. In an attempt to capture its qualitative behavior, however, we have 
devised a toy prescription in terms of a wavelength-dependent diffusion coefficient and measured 
its numerical values for our model galactic disk. In principle, this prescription could be adopted 
as a sub-grid physical process in semi-analytic chemical evolution models as well as cosmological 
simulations. Doing so should help us further constrain the dynamical history of disk galaxies. 
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A. HYPERDIFFUSION AND SHOCK DIFFUSION WITH FIXED MESH 

REYNOLDS NUMBER 

Depending on the system of interest, the Pencil Code requires artificial terms in all dynamical 
equations, except those describing the passive scalar fields, to stabilize the scheme. To simulate 
transonic turbulence with formation of shocks, we include hyper-diffusion and shock diffusion terms 
in our simulations. The hyper-diffusion terms use sixth-order derivatives to damp numerical noise 



at high wavenumber but preserve power on larger scales (Haugen Sz Brandenburg 2004 Johansen 



&z Klahr]|2005 ) , while the shock diffusion terms are of von Neumann type ( Haugen et al.|[2~004 Lyra 



et al. 2008[ ). The usual approach is to set the diffusion coefficients to constant values (^3 and a s 



defined below). However, we have implemented a new strategy to dynamically adjust them so that 
the mesh Reynolds number remains nearly constant. We briefly describe the underlying concept of 
this implementation in this section. 

The hyper-diffusion terms are of the form 

(£f + fa + £f), (A1) 

\ ox b oy b oz b ) 



where Q is the primitive variable to be solved for and 1/3 is the hyper-diffusion coefficient. The 



strength of this operation can in fact be evaluated by comparing Equation ( Al ) with the advection 
term u • VQ. Consider a specific signal (or rather, noise) in Q at wavenumber k. It is damped 
faster than being advected away if 

\u-k\<v 3 (k 6 x + k 6 y + k 6 z ). (A2) 
Since |u • k| < uk < u max /c and k® + ky + k\ ~ k b , where ii max is the maximum magnitude of 



velocity u in the computational domain, Equation (A2) implies ii max ^ v 3 k 5 . We define the mesh 
Reynolds number for hyper-diffusion as 

Re, = ^^, (A3) 

^Nyq 

where k^ yq = it/ max(<5x, 5y, 5z) is the Nyquist wavenumber, and 5x, Sy, and 5z are grid spacing 
in the x, y, and z directions, respectively. The aforementioned criterion for damping signals at 
Nyquist frequency then becomes Re^, < 1. 



Motivated by this criterion, we invert Equation (A3) to find the value of a time-dependent, 
spatially uniform hyper-diffusion coefficient vj, with a fixed mesh Reynolds number Re^: 

^3 = V3{t) = 5 ■ (A4) 
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In other words, we determine the maximum magnitude of the velocity field ii max at the beginning of 
each time step, and use this information to assign for this step the value of the diffusion coefficient v% 



calculated from Equation (A4). This way, we maintain the artificial diffusion at Nyquist frequency 
with roughly the same strength. Due to the high-order dependence of the hyper-diffusion operator 
on wavenumber (~ A; 6 ), the damping of noise is then concentrated at and near the Nyquist frequency 
while quickly diminishing towards longer wavelengths. 

Similarly, we can control the strength of shock diffusion by fixing the appropriately-defined 
corresponding mesh Reynolds number. The shock diffusion terms are of the form V • (V S VQ) 
except the one for the momentum equation, which is written as a bulk viscosity p _1 V (p^ s V • u). 
The diffusion coefficient v s is of the form v s = a s max(— V • u, 0), where a s is a positive constant; v s 
is thus spatially variable and is proportional to the local convergence of the flow. Consider again 
a signal in Q with wavenumber k and compare the strength of shock diffusion with that of the 
advection. One obtains 

|u-k| <v s {k 2 x + k 2 y + k 2 z ) = v s k 2 . (A5) 

Note that this criterion is only meaningful near shock fronts. We therefore define the mesh Reynolds 
number for shock diffusion as 

_ max {(\u x k x \ + \u y k y \ + \u z k z \) /k 2 ) 

-tie s = 7 ~ n j I Ad J 

a s max(— V • uj 

which is the most conservative measurement of the Reynolds number at the strongest local conver- 
gence of the flowQ With this definition, then, we solve Equation (A6) for the value of the constant 



a s with a fixed Reynolds number Re s at the beginning of each time step. 

In all of our simulations, we use Re/j = Re s = 1/4 except for model MTF, in which we use 
Re s = 1/6. 



B. RESOLUTION STUDY 

In this section, we demonstrate that our simulations are numerically converged. The diagnostic 



we choose to present is the y-power -Px-(t, k y ) defined in Equation (22), which is arguably the most 
important measurement from our simulations made in this paper. Figure [12] plots the power 
Pjf it , feinj ) as a function of the advection time i for Aj n j = L and Aj n j = L/8 from model T at 
different resolutions. For the case of Aj n j = L, the curves from resolutions 128x128 to 2048x2048 
are roughly on top of each other, and thus they all exhibit almost the same behavior on the 
two stages of turbulent mixing discussed in Section [4] For the case of Ai n j = L/8, significant 
amounts of power are lost in the low resolution simulations due to their inability to resolve the 



7 The Reynolds number Re s is undefined if there is no position for which V - u < 0, and no shock diffusion operates 
in this case. 
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short wavelength of the injected metal distribution. However, one can see in Figure 12 that the 
difference between each pair of adjacent curves decreases with higher resolutions, and the curves 
for the highest two resolutions, 1024x1024 and 2048x2048, roughly coincide, indicating numerical 
convergence. Therefore, turbulent mixing of metals in our simulations should not be dominated by 
numerical dissipation, and the values listed in Table [3] should be robust. 
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Fig. 12. — Power of the metal tracer field at the injection wavelength Aj n j for Aj n j = L (top) and 
A; n j = L/8 (bottom) as a function of the advection time i in model T. The power is averaged over 
10 snapshots at regular (physical) time interval from t = 6P to 15P. Different lines are the power 
from the same model at different resolutions, and the lines are smoothed by running averages to 
emphasize their general trend with respect to t. 
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Table 1. Adopted Physical Parameters 



Parameter 


Symbol 


Value 


Galactocentric distance 


Ro 


10 kpc 


Angular circular speed 




26 km s _1 kpc" 1 


Orbital period 


p 


240 Myr 


Spiral-arm multiplicity 


m 


2 


Spiral-arm pitch angle 


i 


5.7° 


Inter-arm distance 


L 


3.1 kpc 


Initial Toomre stability parameter 


Qo 


1.5 


Mean molecular weight 




1 


Vertical disk scale height 


H 


95 pc 


Table 2. List of Models 


Model Forcing 3, Equation of State b 


Magnetized Highest Resolution 


Control No Isothermal 


No 


1024x1024 


F Yes Isothermal 


No 


1024x1024 


T No Non-isothermal 


No 


2048x2048 


TF Yes Non-isothermal 


No 


2048x2048 


M No Isothermal 


Yes 


1024x1024 


MF Yes Isothermal 


Yes 


1024x1024 


MT No Non-isothermal 


Yes 


2048x2048 


MTF Yes Non-isothermal 


Yes 


1024x1024 


a If forcing exists, F = 3%; F = 0, otherwise. 




b For isothermal disks, So = 13 Mq pc 


~ 2 and c s ,o 


= 7.0 km s -1 . For non- 


isothermal disks, Sq = 12 Mq pc -2 , c Sj o = 


6.4 km s _1 


, and 7 = 1.8. 



c For magnetized disks, fio = 2. 
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Table 3. Properties of the Mixing Process for Different Metal Tracers in Model T 



Ainj 

(kpc) 


to 
(Myr) 




First Stage 


Second Stage 


(Myr) 


D 

(kpc 2 Gyr- 1 ) 


td 
(Gyr) 


D 

(kpc 2 Gyr" 1 ) 


3.1 


100 


48 


5.2 


0.20 


1.2 


1.6 


41 


18 


3.5 


0.16 


0.38 


0.78 


22 


8.6 


1.8 


0.13 


0.12 


0.39 


12 


4.0 


0.96 


0.11 


0.037 



Note. — A^j is the wavelength of the metal distribution injected from 
the left boundary, to denotes the approximate advection time when the 
mixing process transitions from the first stage to the second, td and D 
respectively represent the decay time constant of the injected distribution 
and the corresponding diffusion coefficient at each stage. 



